Introduction to Astrometry - Exercises

Dora Föhring, University of Hawaii Institute for Astronomy

Aim: To measure the position and uncertainty of a Near Earth Asteroid

0. Prerequisites

If you do not have astroquery installed, you will need to add it to your conda environment (from the command line):

conda install -c astropy astroquery



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
from astropy import wcs
from astroquery.gaia import Gaia
## make matplotlib appear in the notebook rather than in a new window
%matplotlib inline

0.1 Directory Set up & Display Image



In [ ]:

    
datadir = ''
objname  = '2016HO3'



In [ ]:

    
def plotfits(imno):
    img = fits.open(datadir+objname+'_{0:02d}.fits'.format(imno))[0].data

    f = plt.figure(figsize=(10,12))
    im = plt.imshow(img, cmap='hot')
    im = plt.imshow(img[480:560, 460:540], cmap='hot')
    plt.clim(1800, 2800)
    plt.colorbar(im, fraction=0.046, pad=0.04)
    plt.savefig("figure{0}.png".format(imno))
    plt.show()



In [ ]:

    
numb = 1 
plotfits(numb)

1. Centroiding on Images

Write a text file with image centers. Write code to open each image and extract centroid position from previous exercise. Save results in a text file.



In [ ]:

    
centers = np.array([[502,501], [502,501]])
np.savetxt('centers.txt', centers, fmt='%i')
centers = np.loadtxt('centers.txt', dtype='int')



In [ ]:

    
searchr = 5

1.1 Center of Mass



In [ ]:

    
def cent_weight(n):
    """
    Assigns centroid weights
    """
    wghts=np.zeros((n),np.float)
    for i in range(n):
        wghts[i]=float(i-n/2)+0.5
    return wghts

def calc_CoM(psf, weights):
    """
    Finds Center of Mass of image
    """
    cent=np.zeros((2),np.float)
    temp=sum(sum(psf) - min(sum(psf) ))
    cent[1]=sum(( sum(psf) - min(sum(psf)) ) * weights)/temp
    cent[0]=sum(( sum(psf.T) - min(sum(psf.T)) ) *weights)/temp
    return cent



In [ ]:

    
centlist = []
for i, center in enumerate(centers):
    image = fits.open(datadir+objname+'_{0:02d}.fits'.format(i+1))[0].data
    searchbox = image[center[0]-searchr : center[0]+searchr, center[1]-searchr : center[1]+searchr]
    boxlen = len(searchbox)
    weights = cent_weight(boxlen)
    cen_offset = calc_CoM(searchbox, weights)
    centlist.append(center + cen_offset)



In [ ]:

    
print(centlist[0])

2. Identifying Stars in the Field

Ex 1. Write code to identify stars in the field.

One way to do it would be:
Create a new image using an mapping arc sinh that captures the full dynamic range effectively.
Locate lower and upper bounds that should include only stars.
Refine the parameters to optimize the extraction of stars from background.



In [ ]:

    
no = 1
image = fits.open(datadir+objname+'_{0:02d}.fits'.format(no))[0].data

1.a. Create a new image using an mapping arc sinh that captures the full dynamic range effectively. Consider Gaussian smoothing to get rid of inhomogineties in the image.



In [ ]:

    
## Some functions you may want to use
import skimage.exposure as skie
from scipy.ndimage import gaussian_filter



In [ ]:

    
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###

1.b. Create a new image that is scaled between the lower and upper limits for displaying the star map.

Search the arcsinh-stretched original image for local maxima and catalog those brighter than a threshold that is adjusted based on the image.



In [ ]:

    
## Consider using
import skimage.morphology as morph



In [ ]:

    
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###

Plot image with identified stars and target



In [ ]:

    
f = plt.figure(figsize=(10,12))
plt.imshow(opt_img, cmap='hot')
plt.colorbar(fraction=0.046, pad=0.04)
plt.scatter(x2, y2, s=80, facecolors='none', edgecolors='r')
plt.scatter(502.01468185, 501.00082137, s=80, facecolors='none', edgecolors='y' )
plt.show()

3. Converting pixel coordinates to WCS



In [ ]:

    
def load_wcs_from_file(filename, xx, yy):
    # Load the FITS hdulist using astropy.io.fits
    hdulist = fits.open(filename)

    # Parse the WCS keywords in the primary HDU
    w = wcs.WCS(hdulist[0].header)

    # Print out the "name" of the WCS, as defined in the FITS header
    print(w.wcs.name)

    # Print out all of the settings that were parsed from the header
    w.wcs.print_contents()

    # Coordinates of interest.
    # Note we've silently assumed a NAXIS=2 image here
    targcrd = np.array([centlist[0]], np.float_)
    
    starscrd = np.array([xx, yy], np.float_)

    # Convert pixel coordinates to world coordinates
    # The second argument is "origin".
    world = w.wcs_pix2world(starscrd.T, 0)

    return w, world

Find position of Asteroid in WCS



In [ ]:

    
wparams, scoords = load_wcs_from_file(datadir+objname+'_{0:02d}.fits'.format(1), x2, y2)



In [ ]:

    
print(scoords)



In [ ]:

    
wparams, tcoords = load_wcs_from_file(datadir+objname+'_{0:02d}.fits'.format(1), np.array([centlist[0][0]]), np.array([centlist[0][1]]))



In [ ]:

    
print(tcoords)

3. Matching

3.1 Get astrometric catalog



In [ ]:

    
job = Gaia.launch_job_async("SELECT * \
FROM gaiadr1.gaia_source \
WHERE CONTAINS(POINT('ICRS',gaiadr1.gaia_source.ra,gaiadr1.gaia_source.dec),CIRCLE('ICRS', 193.34, 33.86, 0.08))=1;" \
, dump_to_file=True)

print (job)



In [ ]:

    
r = job.get_results()
print (r['source_id'], r['ra'], r['dec'])
print(type(r['ra']))

3.2 Perform Match

Convert Gaia WCS coordinates to pixels



In [ ]:

    
ra  = np.array(r['ra'])
dec = np.array(r['dec'])

xpix, ypix = ### fill in one line here ###

Plot Gaia stars over identified stars in image



In [ ]:

    
f = plt.figure(figsize=(20,22))
plt.imshow(opt_img, cmap='hot')
plt.colorbar(fraction=0.046, pad=0.04)
plt.scatter(x2, y2, s=80, facecolors='none', edgecolors='r')
plt.scatter(xpix, ypix, s=80, facecolors='none', edgecolors='g')
#plt.scatter(xpix[17], ypix[17], s=80, facecolors='none', edgecolors='y')
plt.imshow(opt_img, cmap='hot')
plt.show()

Ex. 2 Find the amount of shift needed.

Match catalogue stars to identified stars and measure amount of shift needed to overlay FoV stars to catalogue.

E.g. Find closest star to one of the Gaia stars near the center of image. Find magnitude of shift. Shift all other Gaia stars and see whether resulting difference is small.



In [ ]:

    
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###
### code here ###

Shift



In [ ]:

    
targshifted = centlist[0] + np.array([xshift, yshift])

Convert shifted coordinate into WCS



In [ ]:

    
wparams, tscoords = load_wcs_from_file(datadir+objname+'_{0:02d}.fits'.format(1), np.array([targshifted[0][0][0]]), np.array([targshifted[0][0][1]]))